
function y=MU_Est_Given_Data_Pref_Mod_Rewt(y, x, cutoffs, b_weights, pi, model, init_theta, p, l_bd, u_bd)

%Given the data, preferences, and model specification, this routine maximizes 
%the forecaster's objective using the simulated annealing (SA) algorithm.

%The sample is assumed to be endogenously startified; the objective
%function is corrected for the assumed unconditional probability that Y=1

%**************************************************************************
%                                   INPUTS
%**************************************************************************

%y:           N x 1 column vector of 1's and -1's;  
%x:           N x ? matrix of covariates
%cutoffs:     N x 1 column vector
%b_weights:   N x 1 column vector
%pi:          scalar; the unconditional probability P(Y=1)
%model:       character string; name of the m-file containing the model m(x,theta)
%init_theta:  row vector; initial value of theta to be used in the simulated annealing algorithm
%p:           scalar; dimension of theta: number of parameters (coefficients) in the model
%l_bd:        (1 x p) row vector of lower bounds for each coefficient
%u_bd:        (1 x p) row vector of upper bounds for each coefficient

%m-files called by this routine: min_sa.m, the m-file designated by the
%variable 'model'

%**************************************************************************
%                                   OUTPUTS
%**************************************************************************

%1 x (p+1) vector: First entry is the value of the maximized empirical
%objective given the data, preferences and model specification. The rest of 
%the vector contains the maximizer (i.e. the maximizing value of theta)


%                           SUBROUTINE STARTS HERE 

% Setting up that part of the objective function that is invariant to the
% model specification (i.e. theta)
%**************************************************************************

N=length(y);
dim=size(x);

%S1 = b_weights.*( y+1-2*cutoffs ); %N x 1

index_y1=find(y==1);
y1=y(index_y1);
x1=x(index_y1, : );
c1=cutoffs(index_y1);
b1=b_weights(index_y1);
N1=length(y1);

S1 = b1.*( y1 + 1 - 2*c1 ); 

index_yneg1=find(y==-1);
yneg1=y(index_yneg1);
xneg1=x(index_yneg1, : );
cneg1=cutoffs(index_yneg1);
bneg1=b_weights(index_yneg1);
Nneg1=length(yneg1);

Sneg1 = bneg1.*( yneg1 + 1 - 2*cneg1 ); 


%Setting the parameters of the SA procedure
%*******************************************

    v=0.5*ones(1,p);
    %initial step size in the p coordinate directions

    c=2*ones(1,p);
    %a 1 x p vector, a controls step size adjustment for each coordinate

    %epsilon=10^(-8);
    epsilon=10^(-6);
    % terminating criterion

    %N_epsilon=6;
    N_epsilon=4;
    % the number of iterations over which the termination criterion
    % must be satisfied for code to stop

    %N_S=25;
    N_S=20;
    % how many searches through the coordinates before adjusting the step size

    %N_T=120;
    N_T=100;
    % how many step adjustments before reducing the temperature
    
    %init_theta=[init_theta; init_theta; l_bd + rand(1, p).*( u_bd - l_bd )];
    dim=size(init_theta);
    n_runs=dim(1);
    % number of runs of the SA algorithm
 
    
    T=[1000  1000  1000];
    r_T=[0.65   0.65   0.65];
    % A number of runs are possible with different starting
    % temperatures and annealing rates (temperature reduction coefficient)
    % in each. In the rth run, T(r) and r_T(r) is used.
          
for r=1:n_runs
   
    mininfo=min_sa( @minus_emp_score_rewt, p, init_theta(r,:), v, T(r),...
                    epsilon, N_epsilon, N_S, c, N_T, r_T(r), l_bd, u_bd );
    %Note: this routine minimizes, so the negative of the objective function is given as an input            

    theta_opt(r,:)=mininfo{1};
    S_opt(1, r)=mininfo{2};

end %SA runs

%Choosing the best result if algorithm ran more than once
S_max=max(-S_opt); 
theta_max=theta_opt(find(S_max==-S_opt), : );
theta_max=theta_max(1,:);

y=[S_max theta_max]

%Below the forecaster's objective is set up as a nested function. Note that
%'minus_emp_score' is specified as function of 'theta' only. This is
%necessary for the SA procedure ('min_sa') to work. The rest of the
%parameters are passed on to 'minus_emp_score' through the outside
%function 'MU_Est_Given_Data_Pref_Mod'.


    function S=minus_emp_score_rewt(theta)

        index1=feval( str2func(model), x1, theta );         %N1 x 1
        index_neg1=feval( str2func(model), xneg1, theta );  %Nneg1 x 1
        
        sign1    =( ( index1-c1 )>0 ) - ( ( index1-c1 )<=0 );
        sign_neg1=( ( index_neg1-cneg1 )>0 ) - ( ( index_neg1-cneg1 )<=0 );
        
        S=(pi/N1)*sum(S1.*sign1) + ( (1-pi)/Nneg1 )*sum(Sneg1.*sign_neg1);
       
        S=-S;

    end %minus_emp_score


 end %MU_Estimation

    
    


